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Abstract. We present a preliminary study of the multipolar structure of 
gravitational radiation from spinning black hole binary mergers. We consider 
three different spinning binary configurations: (1) one "hang-up" run, where the 
black holes have equal masses and large spins initially aligned with the orbital 
angular momentum; (2) seven "spin-flip" runs, where the holes have a mass ratio 
q = M1/M2 = 4, the spins are anti-aligned with the orbital angular momentum, 
and the initial Kerr parameters of the holes j\ = 31 = ji (where j = J/M 2 ) arc 
fine-tuned to produce a Schwarzschild remnant after merger; (3) three "super- 
kick" runs where the mass ratio q = 1, 2, 4 and the spins of the two holes 
are initially located on the orbital plane, pointing in opposite directions. For 
all of these simulations we compute the multipolar energy distribution and the 
Kerr parameter of the final hole. For the hang-up run, we show that including 
leading-order spin-orbit and spin-spin terms in a multipolar decomposition of the 
post-Newtonian waveforms improves agreement with the numerical simulation. 

These are exciting times for gravitational wave (GW) research. Earth-based laser- 
interferometric detectors are collecting data at design sensitivity, and LIGO [1] just 
completed the longest scientific run to date. The space-based interferometer LISA 
is expected to open an observational window at low frequencies (~ 10~ 4 — 10 _1 Hz) 
within the next decade [2]- Following remarkable breakthroughs in the simulation 
of the strongest expected GW sources, the inspiral and coalescence of black hole 
binaries p2 [4j [5] , several groups have now explored various aspects of this problem, 
including spin- precession and spin-flips [6l [7] , comparisons of numerical results with 
post-Newtonian (PN) predictions 02 |10l [TT], multipolar analyses of the emitted 
radiation (HO [12] and the use of numerical waveforms in data analysis p73 | fl4 t fl~5 t 1 1 6J . 

In Ref. [Hj we studied the multipolar distribution of radiation and the final spin 
resulting from the merger of unequal-mass, non-spinning black holes with mass ratios 
q = M1/M2 in the range 1 to 4. The main purpose of this paper is to show preliminary 
results from our attempt to extend the analysis to spinning binaries. 
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A second purpose of this study is to test recent predictions for the spin of 
the black hole resulting from a generic merger. Buonanno, Kidder and Lehner [T7] 
recently introduced a surprisingly accurate model, based on the extrapolation of point- 
particle results, that was shown to be in good agreement with existing numerical 
simulations (see also [H]). An interesting question explored in [17] concerns spin- flip 
configurations. Suppose that initially both black holes have equal Kerr parameters 
(ji = Ji — 32), and spins antialigned with respect to the orbital angular momentum. 
For a given mass ratio q, which value of ji will produce a Schwarzschild black hole? 
These "critical" configurations could be very interesting, since mild variations of 
the parameters around the critical values may produce interesting orbital dynamics 
(eg. spin flips) and complex gravitational waveforms. As argued in [17j . one needs 
unequal masses to be able to produce a Schwarzschild remnant at all. For q = 4 
and zero eccentricity, Ref. [17] predicts that a Schwarzschild black hole should be 
formed when ji = —0.815. A semi-analytical fitting formula [TH] predicts a critical 
spin ji = —0.823. Here we provide a numerical benchmark against which to test these 
analytical models, and possibly other models that may be developed in the future, 
by computing the final Kerr parameter jfj n from a sequence of spinning binaries with 
q — 4 and values of ji G [—0.75, — 0.87] . The numerical results are well fitted by a linear 
relation of the form ja B = — 0.570(j, — 0.842). From this fit, our best estimate for the 
initial spin leading to the formation of a Schwarzschild remnant is ji ~ —0.842 ±0.003. 

The plan of the paper is as follows. In Sec. [T] we introduce the numerical code 
and we list the new simulations considered in this paper. In Sec. [2] we compare 
the multipolar energy distribution of spinning and non-spinning binaries. In Sec. [3] 
we generalize our multipolar decomposition of PN waveforms to include leading-order 
spin contributions in the special case of spins aligned (or anti-aligned) with the orbital 
angular momentum, and we show (in a special case) that the inclusion of spin terms 
improves the agreement with numerical results [19]. In Sec. [4] we study in detail "spin- 
flip" configurations designed to produce a Schwarzschild remnant. Finally, in Sec. [5] 
we discuss the fraction of energy radiated in ringdown in the different simulations. 

1. Numerical simulations 

In this work we compare two sequences of numerical black hole binary simulations. 

Sequence 1 is a series of simulations of non-spinning black hole binaries with mass 
ratio q ranging from 1 to 4. These simulations were performed with the Bam code 
[2"D] , and they were used in [5T] to study gravitational recoil and in [9] to investigate 
the multipolar structure of the emitted gravitational radiation. Sequence 2 consists 
of simulations of binary systems with mass-ratios in the same range, but with non- 
vanishing spins. In particular, we study two families of spinning binaries. For the first 
family, the spin of both holes is either aligned (run uul) or anti-aligned (runs ddl- 
dd7) with the orbital angular momentum. For the second family, the spins of the two 
holes are initially located on the orbital plane, and they point in opposite directions. 
The latter configuration produces surprisingly large recoil velocities [22] , and it is 
sometimes referred to as a "super- kick" configuration (runs ski, sk2, sk4). 

Sequence 2 binaries have been evolved with an advanced version of the Lean code 
described in |23] . The key improvement over the original code is the implementation 
of sixth-order accurate stencils for spatial derivatives, as introduced in [53]. Lean 
is based on the Cactus computational toolkit and uses Carpet [25] for mesh- 
refinement, TwoPunctures [26] for puncture initial data and AHFinderDirect 
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Table 1. Details of the spinning binary simulations of sequence 2. The holes 
start on the x axis with a coordinate separation of 6M for models uul and ski— 4 
and 8M for models ddl-7. Black hole 1 is located at x > and hole 2 at x < 0. 
For the spin-kick (sk) runs Si = — 5*2 and the spins are aligned along the rr-axis, 
while for all other models Si = +S2 and spins are aligned along the 2-axis. A 
common horizon forms at t ca ^ and the number of gravitational wave cycles from 
to = )*cx + 50 M to the peak in the wave amplitude is iV cyc . -Ebmop is the total 
ringdown energy radiated in / = m = 2 and (/ = 2, m = —2); the number in 
parenthesis is the percentage of energy in the (I = 2, m = —2) mode. 
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[27l [28] for horizon rinding. We have used the advanced version of Lean for all 
spinning configurations, except for the model labelled uul in Table [T] which has first 
been reported in Ref. [19]. 

In Table [1] we list the initial physical parameters of all binaries of sequence 2. We 
also list: the time i ca h of formation of a common apparent horizon; the number of wave 
cycles N cyc derived from the phaseQ] of the (I = 2, m = 2) mode; the total radiated 
energy -E^ad and the total radiated angular momentum J ra d, excluding again the initial 
data burst; the energy radiated in ringdown -Eemop; as estimated using the energy- 
maximized orthogonal projection [9], for the dominant modes (I = 2,m = 2) and when 
the symmetry is broken (i.e., for the sk runs) also for (I = 2, m = —2); and the final 
black hole spin jfi n = Jfi n /M| n . For models sk2 and sk4, some angular momentum is 
radiated in the x- and y-directions. This results in a realignment of the final spin, and 
computing jg n becomes more difficult. For this reason, the corresponding entries in the 
table are empty. Following the convention of Ref. [29] , initial data are normalized to 
M = Mi + Mi, while all radiated quantities are normalized to the ADM mass Madm- 
The corresponding details for the non-spinning models can be found in Refs. [9"I l2Tj. 

In the notation of Sec. II E of Ref. [23], the grid setup is {(256, 128, 64, 32, 16, 8) x 
(2,1,0.5), h = 1/80} for the uul run,' {(384, 192, 96, 56, 24, 12) x (3,1.5,0.75), h = 
1/48} for the ski and sk2 runs, and {(384,192,96,56,24) x (6,3,1.5,0.75), h = 1/48} 
for all other simulations. In addition, we have performed higher-resolution simulations 
with h = 1/52 and h = 1/56 of the two models labelled ski and sk2 in Table [1] 

The resulting convergence plot for the (I = 2, m = 2) mode ^4 of the Newman- 
Penrose scalar *5>4 extracted at r ox = 60A/adm is shown in Fig. [1] The differences 
between the high-resolution runs have been rescaled by a factor of 1.72, as expected 

f To remove the initial radiation burst the phase was integrated from t = 50 M + r cx , where r cx is 
the extraction radius, up to the maximum in the wave amplitude. 
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Figure 1. Convergence analysis of the (2 = 2,m = 2) mode of Mj^uM r ^4 f° r 
simulations ski and sk2 of Table [T] The difference between the high-resolution 
runs has been rescaled by 1.72, as expected for sixth-order convergence. 

for sixth-order convergence. 

To estimate uncertainties arising from finite resolution, we use the differences 
obtained for these convergence runs. In order to allow for the fact that the Lean code 
also contains ingredients of lower than sixth-order accuracy, we follow the approach of 
Ref. [S] and apply a Richardson extrapolation assuming second-order accuracy. Using 
these conservative estimates, and the observed l/r cx fall-off of the errors on radiated 
energy and momenta, we obtain very similar uncertainty estimates for simulations 
ski and sk2. At resolution h = 1/48 and extraction radius r cx = 60Madm these 
uncertainties are ~ 4% for the radiated energy and ~ 7% for the radiated angular 
momenta. For the discussion below, it is important to remark that our numerical 
results underestimate the radiated quantities in all cases. 

2. Energy distribution for non-spinning binaries 

The decomposition of gravitational radiation from non-spinning binaries onto spin- 
weighted spherical harmonics (SWSHs) -2^,m can be found in [5], where it was 
obtained by projecting the 2.5PN gravitational waveform^ derived in [3TJ[32]- The 
most notable result of Ref. [5] is that the leading-order term contributing to each 
multipolar component tpi >m of the Weyl scalar ^4 is proportional to the symmetric 
mass ratio r\ = q/(l + q) 2 when m is even, and to r/SM/M — rj(Mi — M<z)/M when 
m is odd. More explicitly: 

5 

e^Mr^i, m = r,Y,9l%(v)(Mn)( g + n y 3 , m even; (1) 

n=0 

§ The SWSH components of h+,h x are related to the corresponding components of ^4 by 
—ih x ) [ = —ipi m /(mQ) 2 , with the exception of 2.5PN contributions to the I = m = 2 
component. It has recently been pointed out that by using the known expressions of the radiative 
multipoles, instead of the waveforms, more information is available and the contribution of each 
multipole can be computed to higher PN order I30I . 



Multipolar analysis of spinning binaries 



5 



71=1 

Here <fi is an orbital phase, including logarithmic corrections in the orbital frequency 
n. The precise definition of this phase, as well as the functional form of the coefficients 

SlmM aiid ^!m( ? ?)> can f° un d in Appendix A of [9]. Each coefficient in the series 
represents a PN contribution of order n/2 to the leading-order (Newtonian) prediction. 




Figure 2. Left: Energy Ei in different multipoles for the unequal-mass binaries 
of sequence 1 as a function of q (from 9 ). Right: Ei for some spinning binary 
configurations belonging to sequence 2, as a function of the multipole index I. 
Continuous (black) lines refer to equal-mass binaries, the dotted (red) line to a 
binary with q = 2, and dashed (blue) lines to binaries with q = 4. 

From Eq. ([1]) it is clear that, for non-spinning binaries, odd-m multipoles are 
suppressed in the equal-mass limit. Since / = m components typically dominate the 
radiation for each given I, this also implies that odd-/ multipoles are suppressed as 
q — > 1. This is confirmed by numerical simulations of the merger. In addition, higher 
multipoles are found to carry a larger fraction of the total energy as q deviates from 
unity. Both of these features are clear from a glance at the left panel of Fig.[2j In [9] we 
also showed that, to leading order, the total energy radiated in the merger scales like 
r] 2 and the Kerr parameter of the final hole scales like 77, providing phenomenological 
fits of these quantities. More general fitting formulas for the final Kerr parameter, 
encompassing also binaries with aligned or anti- aligned spins, can be found in [18] . 

In the right panel of Fig. [2] we show the energy distribution for some of our 
sequence 2 binaries, as a function of the multipole index /. The relative uncertainties 
of the mode energies have been determined in analogy to those of the total radiated 
energy for runs ski and sk2 and strongly depend on the energy content (the "signal 
strength"). We find relative uncertainties of about 5% for values of Ei/Mabm down 
to approximately 5 10 -4 . These grow up to 30% for low signals of 10 -5 ...10~ 4 . The 
figure reveals that even in the presence of spin, odd-l multipoles are suppressed when 
q = 1. As first observed in [33], the hang- up (uu) configuration stays in orbit for a 
longer time and radiates more energy before merging. On the contrary, our spin-flip 
(dd) simulations with q = 4 merge very rapidly (compare the number of cycles iV cyc 
in Table [T|). All dd simulations radiate roughly the same amount of energy, so we only 
show run ddl in the plot. By comparing the ski, sk2 and sk4 runs we confirm that 
our conclusion in [9j, that large-g binaries radiate more energy in higher multipoles, 
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holds true also for these spinning binaries. This is clear from the slopes of E\ as a 
function of I for the different sk runs, and it is nicely illustrated by a comparison of 
the sk2 and sk4 runs. Even though sk4 radiates roughly half the energy radiated by 
sk2, the energy radiated in each I > 2 multipole by the sk4 run is larger. 

3. Leading-order spin-orbit and spin-spin contributions 



For spins aligned or anti-aligned with the orbital angular momentum, the leading-order 
spin contributions to the various multipolar components are most easily derived by 
projecting Eqs. (F24) and (F25) of Ref. [34] onto SWSHs, according to the procedure 
described in [9] . Let Si be the projection of the spin of body i on the axis orthogonal to 
the orbital plane. Si is positive (negative) if the spins are aligned (anti-aligned) with 
the orbital angular momentum. Define the dimensionless spin parameter ji = Si /Mf 
(i = 1, 2) and the spin combinations \ s — (ji+j2)/2, x.a = (ji — j2)/2. Including only 
the dominant spin-orbit (1.5PN order) and spin-spin (2PN) terms, in addition to the 
non-spinning part of Eq. JT]) we find three spin-dependent multipolar contributions: 
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These equations demonstrate that odd-m multipoles do not always vanish for 
equal-mass systems. For example, contains a term which is not proportional to 

the mass difference 5M/M: for equal masses and ji ^ j2, the dominant contribution 
to ^2,1 comes from spin terms. Therefore, by simulating binaries with equal masses 
and unequal spins and by looking at the (I — 2, to = 1) component we can 
disentangle subleading spin effects from the leading-order non-spin contributions, 
and thus facilitate the comparison of spin definitions in PN theory and in numerical 
simulations. Unfortunately, for systems with equal mass and equal spins, such as our 
uul model, f/^i" = 0- However, we can still study the effect of spin terms on the 
convergence of the PN approximation by considering the (I = 2, m = 2) component. 

To see how this is possible, consider first the non-spinning case. By taking the 
modulus of ((T|) we get a PN series relating the orbital frequency and the gravitational 
wave amplitude \ipi m |. The convergence rate of the series to the numerical results can 
be studied as follows. 

First, we observe that the frequency of the gravitational waves lugw m a 
multipolar component ipi >m is related to the orbital frequency f2 by wqw = 
mf2. Therefore, given a time-dependent component ip\ >m of the Weyl scalar, the 
numerical value of the binary's orbital frequency can be estimated as ~ ujom = 
—lm(ipi ym /ipi im )/m [E1IH]- Consider now the modulus of Eq. (fTJ), possibly with the 
addition of spin terms such as Eq. Given |^>j, m (t)| on the left-hand side (as 

obtained from the simulation), at any give PN order we can (at least in principle) 
numerically invert the PN expansion on the right-hand side. Since this expansion is 
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only valid for quasi-circular binaries, in this way we get a post-Newtonian quasi- 
circular (PNQC) estimate of the orbital frequency: ft ~ wpnqc- If t nc PNQC 
approximation works well, wpnqc should be close to 0Jr> m - Furthermore, if the PN 
approximation is converging, the agreement should get better as we increase the PN 
order, i.e. the number of terms in the sum on the right-hand side of Eq. (fTj). 




Figure 3. Convergence of the PNQC expansion. Left: non-spinning binaries; 
right: spinning binaries. In the right panel, thick lines estimate the PNQC 
frequency including the spin terms of Eq. J3), and thin lines omit the spin terms. 

In Fig. [3] we show the relative deviation between wpnqc an d lud2, considering 
the dominant (/ = m = 2) component of the radiation. Let us hrst consider the 
left panel. There we show the relative deviation between wpnqc and lod2 for the 
longest non-spinning, equal-mass simulation considered in [9], at different PN orders. 
At early-times we see oscillations in the relative deviation, that damp away as the 
binary evolves. The magnitude of the relative deviation (wpnqc — ^>D2)l^m\ can 
be taken as an indicator of the accuracy of the PN approximation. The plot shows 
that the convergence of the PN series is not monotonic. The transition from inspiral 
to plunge is roughly marked by the vertical lines, that correspond to the point where 
the orbital frequency equals the innermost stable circular orbit or ISCO (as defined in 
[35] , computed at 2PN and 3PN to bracket uncertainties). At this point the PNQC 
frequency, which only makes sense in the inspiral phase, decouples from lob2- 

The right panel shows the relative deviation between wpnqc & n d u>D2 for our uul 
run. Vertical lines mark again the 2PN and 3PN ISCO, that for such large aligned 
spins corresponds to much higher orbital frequencies: in our case Mfij™Q Spm = 0.247, 
while for zero spins we would get MSIj^q 110 spm = 0.129. Here we also indicate the 
formation of a common apparent horizon (CAH). The fact that the PNQC estimate 
again deviates from u>D2 at the ISCO seems to indicate that the PN notion of orbital 
instability makes sense even for such large values of the spin. The relatively short 
duration of the simulation, and the large wiggles induced by numerical noise, clearly 
illustrate the need to start simulations of spinning binaries at larger separation. It 
is also clear that including spin-terms improves the agreement between the numerics 
and the PNQC approximation at 1.5PN and 2PN orders. The trend is reversed at 
2.5PN, possibly because we are not including 2.5PN spin contributions. Higher-order 
calculations of spin contributions in PN theory and longer simulations will be necessary 
for more accurate comparisons. 



Multipolar analysis of spinning binaries 



8 



4. Producing a Schwarzschild remnant 

A question of particular astrophysical importance concerns the final spin resulting 
from the inspiral and merger of black-hole binaries with arbitrary initial parameters 
jTTJ [TS] . An intriguing special case is that where two black holes with initial spins 
anti-aligned with the orbital angular momentum merge forming a non-spinning hole. 

Our sequence of runs with anti-aligned initial spins and mass ratio q = 4 
(symmetric mass ratio r\ = 0.16) has been designed to bracket the critical point 
of formation of a Schwarzschild black hole, as predicted in Refs. jTTJ [TS]. We 
calculate the final Kerr parameter using energy and angular momentum balance 
arguments, i.e. we compute the final black-hole mass as Mfi n = Madm — -Erad 
and the final angular momentum as Ja n = Jj n i — J ra d- The resulting dimensionless 
Kerr parameter ja n is given in the rightmost column of Table [TJ Applying a linear 
regression analysis to these results and the associated error estimates leads to the 
fitting formula j fila = (-0.570 ± 0.040) [j l - (0.842 ± 0.003)]. A Schwarzschild remnant 
is thus produced when the initial spin has the "critical" value j- m i ~ —0.842 ± 0.003. 
As mentioned above, we generally find the uncertainties due to finite differencing and 
extraction radius to underestimate the radiated angular momentum, so that we expect 
the correct value to be > —0.842. We can compare our result to that of Ref. [18] 
by applying standard error propagation to the uncertainties listed in their Eq. (7). 
Specifically, we solve their Eq. (6) for a, set aa n = and calculate the uncertainty 
Aa 2 = ^ i (9a/9wi) 2 Au 2 , where Vi = S4, S5, t , t 2 and t 3 . The resulting critical 
angular momentum of 0.824 ± 0.019 agrees very well with our value. Both results also 
agree well with that of —0.815 predicted in Ref. [17] . 

Finally, a systematic uncertainty in our numerical results is due to the relatively 
small initial separation of the holes. However, a comparison of non-spinning 
simulations starting at different initial separations shows that most of the angular 
momentum is radiated during the last orbit prior to formation of a common apparent 
horizon [9]. Therefore this systematic error should not significantly affect our results. 
Simulations starting at larger separation are required to verify this expectation, and 
we plan to investigate this effect in future work. 

5. Ringdown energy 

Present ringdown searches in LIGO data are based on matched filtering. For this 
reason, a practical criterion to define the ringdown content of a given waveform is the 
energy-maximized orthogonal projection (EMOP) discussed in [5], which is basically 
matched filtering in white noise. As shown in [9], the EMOP estimate of the energy 
radiated in ringdown is a lower bound on the energy that can be detected by matched 
filtering. Table [1] lists the sum of the ringdown energies radiated in the dominant 
(I = 2,m = 2) and (I = 2,m = —2) components of the radiation. For the sk runs the 
radiation is not symmetric with respect to the z-axis [35], therefore we also list (in 
parentheses) the percentage of £emop radiated in (/ = 2, m — —2). From the data we 
see that as much as ~ 3% of the rest energy of the system is radiated in the ringdown 
phase for equal-mass binaries. From our previous study of non-spinning binaries [9] 
it is reasonable to assume that, for fixed initial Kerr parameters, -Eemop ~ V 2 - Such 
high ringdown efficiencies are good news for the detection of ringdown waves and for 
their use in parameter estimation [371 138] . 

EMOP estimates for runs sk2 and sk4 should be taken with caution. For these 
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runs, the spin axis of the final black hole is tilted with respect to the z-axis of the 
coordinate frame used for the evolutions and for wave extraction, and our chosen 
reference frame is not appropriate to describe the symmetries of the final hole. The 
tilt in the final spin angle produces mode mixing in the ringdown phase: a pure (I, m) 
mode in the frame chosen for numerical computations is a sum of modes with different 
{V , m')'s in the frame whose z-axis coincides with the symmetry axis of the final hole, 
and vice versa. Because of this mode mixing, an estimate of the Kerr parameter of 
the final hole using quasinormal mode fits is not trivial. A detailed treatment of this 
problem will be presented in future work. 
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